Third International Conference on Inverse Design and Optimization in Engineering Sciences 
( ICIDES- III) Editor: C . S .Dulikravich , Washington D.C. October 23-25, 1991. 


N92-13957 


ADJOINT OPERATOR APPROACH TO SHAPE DESIGN 
FOR INTERNAL INCOMPRESSIBLE FLOWS 


H. Cabuk, C.-H. Sung*, and V. Modi, 


Department of Mechanical Engineering, 
Columbia University, ; 

New York, New York 10027 /*H' 


* David Taylor Research Center 
Bethesda, Maryland 20084 

vv 




The problem of determining the profile of a channel or duct (given its upstream cross-section and 
length) that provides the maximum static pressure rise is solved. Incompressible, laminar flow governed 
by the steady-state Navier-Stokes equations is assumed. Recent advances in computational resources and 
algorithms have made it possible to solve the “direct” problem of determining such a flow through a body 
of known geometry. It is possible to obtain a set of “adjoint” equations, the solution to which permits the 
calculation of the direction and relative magnitude of change in the diffuser profile that leads to a higher 
pressure rise. The solution to the adjoint problem can be shown to represent an artificially constructed flow. 
This interpretation provides a means to construct numerical solutions to the adjoint equations that do not 
compromise the fully viscous nature of the problem. This paper addresses the algorithmic and computational 
aspects of solving the adjoint equations. The form of these set of equations is similar but not identical to 
the Navier-Stokes equations. In particular some issues related to boundary conditions and stability are 
discussed. The use of numerical solvers is validated by solving the problem of optimum design of a plane 
diffuser. The direct as well as the adjoint set of partial differential equations are discretized using a finite- 
volume formulation. Each of the resulting set of algebraic equations are then solved numerically to obtain 
a change in profile that will ensure an increase in the static pressure rise. Upon successive applications 
of this procedure, an “optimum” profile is obtained beginning with an initial guess of a diffuser profile. 
Such optimum diffuser profiles are obtained at Reynolds numbers varying from 10 to 2000. The optimality 
condition, that the shear stress all along the wall must vanish for the optimum diffuser, is also recovered 
from the analysis. It is shown that numerical solutions obtained in this fashion do satisfy the optimality 
condition. 

1. INTRODUCTION 

A shape optimization problem is one in which an objective function defined on a domain and/or on its 
boundary through the solution of a boundary value problem, is minimized (or maximized) with respect to 
the variation of the domain. One problem of this nature is “What is the shape of a body (of given volume) 
which has minimum drag when moved at constant speed in a viscous fluid?”. Pironneau (1973) addressed 
this problem in Stokes flow for a three-dimensional unit-volume body. It was shown that at optimality the 
normal derivative of the velocity is constant along the boundary of the body. In addition it was also shown 
that the general shape of the body is similar to a prolate spheroid including a conical front end and rear ends 
of angle 120 degrees. However, due to the lack of a numerical Stokes flow solver, a complete body profile 
could not be obtained. 

In a subsequent study, Pironneau (1974) derived the change in energy dissipation due to a small hump 
on a body in uniform, steady, laminar flow. Using the above result in conjunction with variational methods 
of optimal control “necessary optimality conditions” for four minimum-drag problems were obtained. These 
conditions lead to a set of equations for an additional set of variables called the “co-state” or the “adjoint” 
variables as opposed to the “direct” variables which are the unknown velocities. At the time Pironneau (1974) 
was unable to carry out such a numerical integration. Instead, however using a boundary layer assumption 
he was able to prove that a two-dimensional unit-area body with the smallest drag has a wedge-shaped 
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front end. In a subsequent work Glowinski and Pironneau (1975) presented numerical computations of the 
minimum-drag profile of a two-dimensional body in laminar flow, although with a Reynolds number large 
enough (between 1,000 to 100,000) to permit a boundary layer approximation. The present study belongs 
to this class in its theoretical approach with particular emphasis on computation of optimum profiles in the 
absence of simplifying assumptions such as Stokes flow or thin boundary layers. 

Another related class of optimum design problems is the question of determining the profile of a two- 
dimensional body that will attain a desired surface pressure distribution. The body is assumed to be in 
otherwise uniform flow. The designer usually has a better understanding of how the performance is related 
to the the pressure distribution than the relationship between the profile and the performance. In recent 
survey paper, Jameson (1988) suggests that the design problem be treated as a control problem in which 
the control is the profile of the boundary. He also provides a comprehensive summary of the earlier related 
studies in this direction. In a significant step towards addressing real flows Giles et. al. (1985) addressed 
the problem of shape design for flows governed by the two-dimensional Euler equations. They write the 
two-dimensional Euler equations in a streamline coordinate system and for fixed pressure distribution obtain 
a Newton solution for the unknown surface coordinates. 

In the present study optimum design of an internal flow component such as a diffuser in laminar flow 
is considered. The problem of determining the profile of a plane diffuser (of say, given upstream width and 
length) that provides the maximum static pressure rise is formulated using a variational method derived 
from optimal control theory. Careful consideration of the numerical stability of the adjoint equations we 
have been able to demonstrate the feasibility of optimum design in the context of laminar Navier-Stokes 
equations without the additional boundary layer assumption. 


2. STATEMENT OF THE PROBLEM 

Consider a plane diffuser as shown in figure 1 of given upstream width W x and given length L with 
incompressible, laminar flow through it. The flow is governed by the incompressible, steady forms of the 
Navier-Stokes and continuity equations. These are: 


u j u t\y = + Mijj v 7 

where p* = p/p. Here tq, p, p, and u are the velocity components, pressure, density and kinematic viscosity 
respectively. 

A no slip condition is imposed on the bounding wall. Dirichlet type boundary conditions are assumed 
at the entrance and exit, specifically, it is assumed that the streamwise velocity component at the entrance 
and exit is specified and the transverse velocity component at the entrance and exit is assumed to be zero. 
Symmetry conditions are assumed at the centerline. All velocities and lengths are scaled using the average 
entrance velocity, V, and the diffuser entrance width W x throughout the paper. Hence the Reynolds number 
for the flow through the diffuser is defined as Rc = (V • W x ) ju . 

It is desired that the optimum diffuser profile be such as to maximize the value of this parameter for a 
given upstream width and length. Since pressure may vary across the diffuser inlet and exit regions it was 
decided to choose the change in the flow weighted integral (over the exit and inlet cross-sectional areas) of 
the static pressure rise as the objective function. This quantity is given by: 


J(Tm) = / p*u,-n,*<fs + / p'uitiids (2) 

Jr, Jr 0 

where n* is the t th component of the unit normal vector and Tm is the portion of the diffuser wall that is 
to be shaped. The goal then is to determine the diffuser profile that maximizes the above function. The 
normalized diffuser length, L/W u (henceforth simply called the length) is kept constant. The normalized 
exit width W 2 /W x ^ (henceforth simply called the exit width) is left arbitrary, and its actual value for the 
optimum diffuser is part of the solution to the problem and is determined along with the rest of the profile. 
Since the only mechanism for total pressure drop in the diffuser is viscous dissipation, the optimum profile 
is also the profile for which the viscous dissipation is a minimum. 
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3. MATHEMATICAL FORMULATION 

In this section, the variation of the objective function with respect to the variation of the boundary is 
obtained by means of a perturbation type of analysis. This analysis follows from arguments not unlike those 
used for optimum design in potential flow, in an earlier paper by (Jabuk and Modi (1990). 

First the variation of the solution of the direct problem due to boundary variation is obtained. Let p(s) 
be an arbitrary function of arclength s, defined on Tay, and let € be a positive number. Here Tw is part of 
the boundary that is to be shaped. The whole boundary, including the wall of the diffuser, the centerline 
and the inlet and exit areas, is denoted by T and the domain enclosed by T is denoted by Q. Let each point 
on Tjvf be moved by ep{s) along the outer normal direction. The curve constructed in this way is denoted 
by r M , € and the new domain is denoted by fl e as shown in figure 1. Let (uj,p e ) be the solution of (1) in the 
new domain Let (4> X) x) be defined as follows. 


4>i = lim € 1 [«? — u t ] € ft , 

f — -0 

n = lim e~ 1 [p* — p*l 6 0 . 
«— *0 


( 3 ) 


Then (u^p*) can be written as: 


u* = u, + e4>i 
p e = p* + ex 


( 4 ) 


Since both (u*,p') and (u t ,p*) satisfy the Navier-Stokes equations, it can be shown that (<f>i,x) satisfy the 
following set of equations: 


Ki = o 

u j4>tj + = — ir.t 9- i/<f>i jj 


( 5 ) 


In a similar way it can be shown that on the fixed portions of the boundary 


4> { =0 on (T — T m ) 


(6) 


since both u* and u t satisfy the same boundary conditions. 

The next step is to derive the conditions satisfied by on Tm. Consider a point P on Tat, and a 
corresponding point P t on Ta/.c such that P t lies on the outward normal n, as shown in figure 1. Assume 
that ep(s) is positive. A Taylor’s series expansion of u* about the point P , evaluated at £ = ±\ along the 
normal direction n is 

= < 1 P + {^) p + °^ 2 ) 

= u*] F 4- t4>i\p + 

Since the velocities satisfy the no slip condition on (i.e. u * ] = u x ] p = 0); 




on r w . 


( 8 ) 


The first variation of the objective function is obtained next. The value of the objective function for the 
new domain is given by, 

J(r Ml .) = [ p‘u‘nids+ [ p*rx‘nida (9) 

Jr, Jr 0 

The first variation of the objective function, 6J f is defined by the relation 


J(r M ) = eSJ + 0(e 2 ). 


( 10 ) 


The first variation of the objective function can be shown to be: 
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which is an integral expression over the entrance and exit boundaries. The next step is the transformation 
of this integral from one that is over T/ and T 0 to one that is over This is achieved through the 

introduction of an adjoint variable problem. The inner product of the perturbation equations ( 5 ) and the 
adjoint variables, (zi,r), integrated over the domain, and added to (ll), and after using the divergence 
theorem gives 


^ J = ^ 7T (Ui - Zi) Kids + V j> ~ ^ 9 

+ j> {r<f>i :n, - &z,*uyny - ds + J J nz x ^dA 

j J <j>\ + t ijZ x j ■+■ UjZj'i r tl ) dA . 


( 12 ) 


The adjoint problem has to be defined such that the domain integrals in (15) vanish identically. The choice 
of boundary conditions for these equations is made such that the only nonzero terms are those that are 
integrals over r^, the wall that is to be shaped. Let us define the following adjoint problem 


Zi y i = 0 in ft 

vzijj -I- uy(*t,y + fj,i) - r , = 0 in n 

Zi = u, on T. 


Using ( 6 ), ( 8 ), and (13), equation ( 12 ) can be written as 


SJ 


!*>(£)(£)*• 


(13) 


(14) 


In the above equation, the integration is over the boundary that is to be shaped. We can choose p(s) as: 


p(s) = u(s) 



(15) 


since that would ensure a positive change in the objective function, J, for a sufficiently small non-negative 
weighting function, a>(s). The function p(s) provides the boundary movement for a positive change in J. 
To evaluate p(s) we need to solve the direct problem (i.e. Navier-Stokes equations) given by (l) and the 
adjoint problem in z, given by (13). Note that the optimality condition is satisfied when either the shear 
stress, duijdn , or the adjoint shear stress, dz x /dn ) on the walls vanishes. The former criterion for optimum 
diffuser profiles was also pointed out by Chang(l976). 

It will be shown that the above formulation is equivalent to the earlier work of Glowinski and Pironneau 
(1975). By a change of variable, the adjoint problem can be transformed into the following form: 


=0 in 0 

VWijj + uyt + u>/uy fi - - <?,i = -vyu,' t y in ft 

= 0 on r. 


(16) 


where 2w x — (zi — u,) and 2 <7 — — p* + (l/ 2 )u^ — 2 uyWyJ. The first variation of the objective function 

then becomes 



The form of the adjoint variable problem defined by (16) is identical to that derived by Glowinski and 
Pironneau (1975). Either one of the above adjoint problems can be solved numerically to obtain the next 
shape. However upon examination of (16), it becomes evident that the u*yuy,« term may lead to a numerically 
unstable scheme. This is because the approach to steady state would be attained via an iterative “time 
evolution” like scheme that would then be of the form dw/dt = u;(const) + • • This form is likely to result 
in the exponential growth of the inevitable roundoff and truncation errors present at any iterative step. Also 
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the presence of the inhomogeneous term, -u ; u, ,, in the above equations may lead to a linear growth of the 
roundoff and truncation errors in the numerical computations. It is expected that these numerical difficulties 
will be absent in the (z { ,r) formulation of the adjoint variable problem obtained in this paper and given by 
(13). Hence this is the set of equations for which the algorithm for the numerical solution of the adjoint 

problem is developed. , , 

As pointed out by Pironneau (1974), the adjoint equations do not seem to arise from any identifiable 
physical phenomenon. It is however possible to demonstrate that the adjoint variable problem is associated 
with a certain artificially constructed flow. A change of variables leads to the following form. 

0 in f) 

0 in n ( 18 ) 

uj = — u,- on T. 



where z\ 
Compare 
form. 


= -z,, u' = -u,, and r' = -r. The first equation in (18) is identical to the continuity equation, 
the second equation in (18) with the Navier-Stokes equation written here in a slightly different 


- Uj (u t y + Uy ti ) -p.i = 0 


where p = p* - (l/2)u*. Observe that the problem in adjoint variable z\ is analogous to the the Navier- 
Stokes problem in variable u, with the following exception: the convective velocities in the adjoint problem 
are specified rendering the problem linear and are obtained from the direct problem. These convective 
velocities, u', are identical in magnitude but opposite in direction to those of the “direct” problem. The 
boundary conditions for the adjoint variables are z\ = -u, on T. Hence on the walls they imply a no 
slip condition as in the direct problem. But at the inflow and outflow boundaries, “adjoint” flow ui found 
entering at the domain exit To and leaving at the domain entrance T/, thus suggesting an “adjoint” flow in 

the direction opposite to that of the actual flow. ^ 

The above interpretation of the adjoint variable problem will be useful in constructing a modified 
problem whose solution will provide numerical values, albeit approximate, for the shear stress, du i/dn, and 
the adjoint shear stress, dzjdn , in (15). It is found that a shape optimization algorithm that obtains its 
boundary movement from these approximate numerical solutions does indeed lead to diffuser shapes that 
satisfy the optimality condition. 


4. NUMERICAL ASPECTS 

The boundary conditions chosen for the diffuser in the above formulation are of Dirichlet type. A 
parallel flow of arbitrary distribution is assumed to exist at the diffuser entrance and exit. These boundary 
conditions are clearly unrealistic both from a practical as well as computational standpoint. However this 
is the only set of boundary conditions for which we have been able to derive the adjomt variable problem. 
Given this limitation it was decided to verify whether the boundary movement suggested by (15) would 
continue to provide a means to obtain optimum shapes even if some of the Dirichlet conditions were replaced 

with computationally suitable Neumann conditions. 

Boundary Conditions for Navier-Stokes Equations 

A parallel flow assumption at the upstream boundary implies Dirichlet boundary conditions for both 
the velocity components. Instead a computationally desirable Neumann condition for the transverse velocity 
component (du 2 /dn = 0 on T;) is substituted while retaining a Dirichlet condition for the streamwise com- 
ponent A parabolic profile corresponding to a fully developed laminar flow is specified for this component. 
At the downstream boundary the parallel flow assumption is replaced with computationally desirable Neu- 
mann conditions for both the velocity components ( du^dn = 3u 2 /3n = 0 on To). Similar approximations 
will be made in the solution of the adjoint variable problem, keeping in mind the reversal of the role of 
entrance and exit boundaries. At the solid wall, a boundary whose profile is to be determined, a no slip 
condition is enforced. At the diffuser centerline the usual symmetry conditions are used since the flow is 
assumed to be symmetric. At the entrance, exit, and wall, pressure has been extrapolated from the within 
the domain by assuming that the second derivative of the pressure vanishes on the domain boundary. At 
the centerline symmetry condition is imposed for the pressure. 
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Boundary Conditions for Adjoint Equations 

The role of entrance and exit are reversed for the adjoint equations. Therefore, at the exit boundary 
a Dirichlet type condition is used only for the streamwise component of the co-state vector. Therefore we 
set Z\ — tii on To, with ui taken from the solution of the Navier-Stokes equations. For the remaining 
component z 2 of the co-state vector at exit and for both components of the co-state vector at entrance, 
Neumann conditions are employed instead. At the wall where all velocity components vanish and therefore, 
Zi , the co-state vector that is analogous to the velocity is set to zero. The adjoint variable, r* , is analogous 
to the pressure term in the Navier-Stokes equations and hence no analytical boundary condition is available 
for this variable. However, a computational boundary condition is implemented for this variable. The value 
of r* is extrapolated to the boundary from values at interior points assuming that the streamwise second 
derivative vanishes at the boundary. This is done at all boundaries except at the centerline where asymmetry 
condition is enforced. 


5. NUMERICAL SOLVERS 

Navier-Stokes Equations Solver 

The primitive variable form of the incompressible steady Navier-Stokes equations is solved using an 
artificial compressibility formulation due to Chorin(l967). In this formulation, the continuity equation 
is modified using the time derivative of the pressure term. The steady-state solution of the Navier-Stokes 
equations is then obtained as the large time solution of the unsteady momentum equations with the perturbed 
divergence equation. These unsteady equations are: 


Vt + fPui'i = 0 

u, it + (uyu.) } . = -p* + IsUijj 


(20) 


where f) is analogous to the speed of sound. Note that these equations do not represent any transient physical 
phenomenon and hence the transient solution has no physical meaning until steady state is attained. This 
is indicated by the vanishing of the time derivative terms in the numerical solution. 

The equations are normalized using the velocity and length scales V and W\ defined earlier. In addition 
time and pressure are normalized using the ratio W x /V and pV 2 respectively. The Reynolds number of the 
flow through the diffuser is then given by Re — (V W x ) / v. 

The equations are discretized in space using a finite volume formulation. The spatial discretization is 
performed on the conservative form of the governing equations using a central difference scheme. 

An explicit one-step multistage Runge-Kutta stepping scheme is used for integration in time. Since 
transient behavior is not an issue and a larger time step is desirable, a four-stage Runge-Kutta scheme 
with first order accuracy in time and a relatively high Courant-FViedrichs-Lewy (CFL) number has been 
chosen. In order to improve the convergence rate, a local time step is computed for each cell at each elapsed 
time level. These time steps have been estimated from a stability analysis of the algorithm. A fourth 
order linear artificial dissipation term is introduced to damp the high-frequency oscillations associated with 
the so-called sawtooth or plus-minus waves, i.e. waves associated with the shortest wavelengths. Implicit 
residual smoothing is performed at each iteration to enhance the stability region of the technique. A more 
complete discussion of the finite volume formulation, stability considerations, local time stepping, artificial 
dissipation, implicit residual smoothing and the computational boundary conditions is provided in Cabuk, 
Sung and Modi (1991). 

The computational grid is generated by solving a set of elliptic partial differential equations similar 
to those suggested by Thompson et al. (1974). The set of algebraic equations thus obtained is solved by 
successive over-relaxation (SOR). A typical grid is shown in figure 2. Grids generated by this method were 
nearly orthogonal and the cell dimensions in each direction are approximately equal. 


Adjoint Equation Solver 

The solution to the adjoint set of equations is obtained as the steady state solution to the following set 
of equations: 


r t = -/9 2 2i,. 


z;,< = vzi.jj + tiy (*,y + *y.<) - - (z/tz*),, - r* 


( 21 ) 
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where r* = r — A nonlinear term ^ (zk^k ) , is introduced in the above equation to enhance the rate 

of convergence. The utility of this term was established by means of preliminary calculations performed on 
a straight duct geometry where an exact solution of the Navier-Stokes solution is known for fully developed 
laminar flow. 

The equations are normalized following a procedure similar to that utilized for the Navier-Stokes equa- 
tions. The nondimensional form of the above equation is identical to the equations above with the exception 
of the first term on the right hand side of (21) where the kinematic viscosity, i/ 1 is replaced by the reciprocal 
of the Reynolds number. 

The numerical algorithm for the solution of the adjoint set of equations is essentially similar to the 
algorithm for the Navier-Stokes equations. Some subtle but important differences do exist since the equations 
solved are after all not the same. A discussion of the numerical algorithm is presented here, since this solution 
to the best of our knowledge represents the first successful numerical solution of the adjoint set of equations 
in the absence of either a thin boundary layer or a Stokes flow assumption. Spatial discretization is carried 
out by centered-difference finite volume formulation. The term, uy (z^y + Zji ) , on the right hand side of (21) 
is not in a divergence form. In the treatment of this term the velocities, uy, which have already been obtained 
by the Navier-Stokes solver, are treated as known quantities and are assumed constant inside each cell. Hence 
the volume integral over the cell is performed by applying the divergence theorem to the remaining part of 
this term, i.e. (z l( y + zy,»). 

The other terms in (21) are treated in the same fashion as the finite volume formulation of the Navier- 
Stokes equations. A fourth order linear artificial dissipation term is introduced to damp high-frequency 
oscillations. Time integration is carried out by a Runge-Kutta scheme with local time stepping. The 
discrete form of the equations for the adjoint problem are: 

- (ASj + BSj + CS K )q = E^-(SI6f + SJSj + SK6 7 K )q - eK(6 } + 6j + S^q (22) 
at AK 


where 


and 


with 


’ 0 

p 2 six 

fi 7 SIY 

p 2 siz 


SIX 

£/ + (u! - Zl)SIX 

(u x -z x )SIY 

(ui — z 1 ) SIZ 


SIY 

(u 2 - z 3 ) SIX 

U+(u 2 -z 2 )SIY 

(u 2 — z 2 ) SIZ 


. SIZ 

{U3- 23) SIX 

(u 3 -z 3 )SIY 

U + (u 3 - Z 3 ) SIZ. 


* 0 

P 2 SJX 

p 7 SJY 

P 7 SJZ 


SJX 

K+ ( Ul - Zl )SJX 

(u x - zi)SJY 

(ui — zi) SJZ 


SJY 

(u 2 — z 3 ) SJX 

V +[u 2 -z 2 )SJY 

(u 2 — z 2 ) SJZ 


. SJZ 

(u 3 - z 3 ) SJX 

(u 3 - 23) SJY 

V + (u 3 - z 3 ) SJZ. 


0 

P 7 SKX 

P 7 SKY 

p 7 SKZ 

1 

SKX 

W + (u! - z^SKX 

(uj - Zl )SKY 

(ui — zj) SKZ 

SKY 

(u 2 -z 3 )SKX 

W + (u 2 - z 2 )SKY 

(u 2 — z 2 ) SKZ 

SKZ 

(u 3 — z 3 ) SKX 

(U3-Z 3 )SKY 

W + (u 3 - z z )SKZ. 


0 0 0 01 


0 0 10 
Lo 0 0 1J 

9 = [r*,z i,Z 2 ,z 3 ] T 

U = mSIX + u 2 SIY + u 3 SIZ, SI = SIX 7 + SIY 7 + SIZ 7 , 

V = UiSJX + u 2 SJY + u 3 SJZ, SJ - SJX 7 + SJY 7 + SJZ 7 , 

W = uiSKX + u 2 SKY + u 3 SKZ, SK = SKX 7 + SKY 7 + SKZ 7 
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The volume of the cell is AV and (SIX, SIY, SIZ), (S J X, S JY, S J Z) and (SKX, SKY, SKZ) are the 
surface-area vectors normal to the /, J and K cell surfaces, respectively. SI, SJ and SK are the squares 
of the surface areas of I, J and K cell surfaces, respectively. The first, second and fourth order centered 
differences are defined in the same fashion as in Sung(1987). The maximum local time step permitted for 
stability is obtained by neglecting both the viscous and the artificial dissipation terms in the adjoint problem 
and is given by 

At < CFL J (23) 

The maximum eigenvalue, Aq, in the above equation is estimated as 


A 



U + y/O 2 + 0 2 C 2 


(24) 


where 

V=\V\+ I SIX (u x -z x )| + I SIY (u 2 - z 2 ) | + I SIZ (u 3 - 23) | 

+ |V| + \SJX (u x — z x ) | + |SJK (u 2 — z 2 ) | + \SJZ (u 3 — 23 ) | 

+ \W\ + |SATX (u, - zi) | + | SKY (u 2 - z 2 ) | + \SKZ (u 3 - z 3 ) | 

and 

c 2 = (|5/x| + |s/r| + \siz\) 2 

+ (|5JX| + \SJY\ + \SJZ\) 2 
+(\SKX\ + \SKY\ + \SKZ\) 2 . 

The maximum eigenvalue of the resulting matrix system, including both the viscous terms and the artificial 
dissipation term has been estimated as 


A 0 =y/xj + (4fie- 1 5//AV + 16 ctf ) 2 

+y/X i J + {4Re~ l SJ/AV + 16 7k) 2 (25) 

+y/x 2 K + (4ffe~ 1 Sif/AV r + 16 7k) 2 


where 


A J 
A j 

A k 


1 

2 

1 

2 

1 

2 


U K + y/&%T^SK 


and 

Oi =\U\ + I SIX (u x - z x ) | + \SIY («2 -Z 2 ) I + I SIZ(u 3 - z 3 ) I 
Uj =\V\ + \S JX (u x - z x ) | + |SJK (u 2 - z 2 ) I + I SJZ (u 3 - z 3 ) | 
tJ K =|W| + | SKX (u x - z x ) | + | SKY (u 2 - z 2 ) | + | SKZ (u 3 - z 3 ) | 

Then the local time step has been computed from (23) with the maximum eigenvalue given by (25). 


Profile Modification Algorithm 
The principal steps of the optimization procedure are; 

a) Choose an initial diffuser profile. 

b) Generate a computational grid that conforms to the diffuser wall. 

c) Obtain the steady state solution to the direct problem. 
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d) Obtain the steady state solution to the adjoint problem, by treating the required velocities as known 
from step (c). 

e) Compute du x /dn and dzi/dn from the solutions in steps (c) and (d) respectively. Choose a non-negative 
weighting function u>(s) and hence obtain p(s) from (15). 

f) Move nodes on the diffuser wall to be profiled along the outer normal direction by p(s). The curve 
connecting the nodes after this movement represents the new diffuser profile. 

g) Go to step (b) unless the change in diffuser pressure rise obtained from step (c) is smaller than a desired 
convergence parameter. 

The iterative profile modification process is continued until the change in pressure rise is a small fraction 
of the total pressure rise. An alternate method is to continue the process until the value of p(s) everywhere 
along the wall is less than a critical value. In step (e), the weighting function, u;(s), is chosen to be 
proportional to the arclength, s, along the diffuser wall measured from the diffuser entrance. This ensures 
that the entrance width is maintained constant but the exit width may vary with the diffuser profile. When 
shifting the diffuser wall to a new curve obtained from step (f) some care must be exercised since the curve 
is being redefined using only a finite number of discretely spaced points. Checks are performed on the 
location of points on the new curve to ensure that boundary nodes do not conglomerate or coalesce after 
their movement to a new position. Heuristic measures are also adopted to ensure that the appearance of 
small amplitude wiggles in the new profile are damped to some extent so as to prevent the growth and built 
up of numerical errors in the subsequent calculation 

RESULTS AND DISCUSSION 

Using the numerical solvers and the profile modification algorithm described above, optimum diffuser 
profiles have been obtained for a single diffuser length L/W l = 3 at Reynolds numbers Re=50, 100, 200 and 
500. A sound speed, /? 2 , of 2 for the Navier-Stokes equations and 2.5 for the adjoint equations was used at 
all Reynolds numbers. The calculation at Re=200 (henceforth called the reference case) has been examined 
in particular detail to establish issues of convergence and accuracy. 

The reference ca se was first examined for convergence of the profile modification algorithm. For this 
purpose, a computational grid of 61 nodes in the i and 31 nodes in the y directions is employed, both for 
the Navier-Stokes as well as the adjoint variable problem. Beginning with an initial shape the diffuser profile 
was obtained after each application of the shape modification algorithm. The initial profile and some of the 
intermediate profiles are shown in figure 3. The change in the profile shape is observed to be small between 
the fourth and the ninth iteration and the change was found to be insignificant after nine iterations. Hence 
the iterative process is stopped at the ninth iteration providing a reasonably converged optimum shape. The 
question of computational accuracy of the solvers and hence the accuracy of the optimum profile is addressed 

next. . . , 

The precise error due to a finite grid size on the optimum profile is difficult to determine since the actual 

optimum curve is not known apriori, nor are any other calculations or experimental data available. However 
one way to estimate the effects of the unavoidable truncation errors in a numerical calculation is to obtain 
the optimum diffuser profile using progressively finer grids until the there is no change with grid size. Once 
again the reference case of Re=200 was examined for this purpose using grids of 31 by 16, 61 by 31 and 
finally 121 by 61. The optimum profiles obtained using the three grids are shown in figure 4. The results 
show that the difference between the shapes is negligibly small, providing some evidence that at these grids 
the contribution of the truncation errors may not be significant. In view of this observation, a grid size of 
61 by 31 is found to be a suitable compromise between accuracy and computational work for the results 
presented here. 

In an earlier section we proposed that it was computationally desirable to replace some of the Dinchlet- 
type boundary conditions with Neumann-type conditions in both the Navier-Stokes and the adjoint equation 
solvers. To justify at least partially the validity of solving the modified numerical problem we must verify 
whether the optimum shapes obtained in this fashion do indeed satisfy the optimality condition, i.e. vanishing 
shear stress on the wall, arising from the analysis. 

In figure 5 the wall shear stress normalized by the corresponding value for a straight duct, is shown 
for the optimum shape as well as at several intermediate stages of iteration. The wall shear stress for the 
optimum shape is found to be vanishingly small for all but 10 percent of the wall at the upstream end. 
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The shear stress distributions at intermediate iterations demonstrate a monotonous decrease towards the 
optimum values. Closer examination of the shear stress for the optimum and intermediate shapes at other 
Reynolds numbers confirm the same behavior as well. Hence the results obtained do provide some aposteriori 
justification for the boundary condition approximations made in the modified numerical problem. 

Further justification is sought by examining the behavior of the objective function for the reference 
diffuser. The velocity averaged static pressure rise (i.e. the objective function defined by equation 2) is 
shown in figure 6 at successive iterations of the shape modification process. The objective function for this 
modified numerical problem is indeed found to increase with each application of the boundary movement 
suggested by equation (15). The area averaged static pressure rise through the reference diffuser also increases 
with with shape modification as seen in figure 6. These observations are found to be valid at calculations at 
other Reynolds numbers in the present study as well. 

In addition to the reference case, calculation of the optimum diffuser profile was carried out at three 
other Reynolds numbers, Re= 50, 100 and 500. In figure 7, these profiles are shown for a diffuser of L/Wi = 3 
for a grid of 61 by 31. At lower Reynolds numbers the optimum diffuser profile permits a larger exit area to 
inlet area ratio as one would expect higher viscous effects to support greater diffusion without separation. 
The angle at which the diffuser profile departs at the upstream corner is difficult to compute accurately since 
the flow in that corner may not be accurately resolved. Nevertheless, the approximate angle decreases from 
56 degrees to 19 degrees as the Reynolds number increases from 50 to 500. For the Reynolds number range 
in which numerical solutions are presented here, further refinement of the grid did not lead to any significant 
change in the optimum profile. This was not found to be true of computations at Reynolds numbers higher 
than 500. 

To evaluate the performance of the optimum diffuser, a pressure recovery coefficient, C p , is defined, 
which is the ratio of the static pressure rise of the optimum diffuser to the static pressure rise for an ideal 
diffuser (in potential flow) with the same W 2 /IV 1 ratio as the optimum diffuser. Note that the denominator 
of this ratio is independent of the actual profile between the upstream and downstream cross-sections of the 
diffuser. Using C p as a parameter, the performance of the optimum diffuser is now compared with that of 
a straight walled diffuser with the same W 2 /W 1 ratio at several different Reynolds numbers in the laminar 
regime. The C p values of straight diffusers are found numerically using the Navier-Stokes solver on the 
straight walled geometry without any shape modification steps. As seen from figure 8, the C p values for the 
optimum diffusers are always higher than those for straight diffusers. 

In conclusion, the feasibility of shape optimization for incompressible laminar flows has been demon- 
strated. This approach may also be adopted to other domain optimization problems where the performance 
depends on the geometry of the component, and flow is governed by the viscous laminar flow (either com- 
pressible or incompressible) equations. It may also be possible to consider variations of the objective functions 
depending upon the design criterion of interest. All computational results presented in this paper were car- 
ried out either on an Intel 30386 33MHz microprocessor based machine or on a micro VAX II workstation. 
The CPU times for these calculations are of the order of several hours. The theoretical framework as well as 
numerical solution code for the extension of the method to three-dimensional flow now exist and such flows 
are the subject of study by the authors at present. 

The research reported in this paper is based upon work supported by the National Science Foundation 
under Grant No. CBT-87-10561. 
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Figure (1) Schematic diagram of a plane diffuser. Flow enters at upstream boundary T/ and 
exits at the downstream boundary r 0 . The wall to be shaped is T M and symmetry 
line is Tc- 



Figure (2) A typical computational grid for a plane diffuser obtained using the grid generation 
program. Grid size is 61 by 31. This was the domain for the optimum diffuser at 
Re=200 and L/Wy = 3. 
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Figure (3) Profiles of a reference diffuser at successive iterations. The grid size is 61 by 31. 0 
: Initial shape, □ : First iteration, A : Fourth iteration, * : Ninth iteration. 



Figure (4) Effect of grid size on optimum profile of a reference diffuser. O : 31 by 16, n : 61 
by 31, * : 121 by 61. 
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Figure (5) Normalized wall shear stress at successive iterations for a reference diffuser. The 
grid size is 61 by 31. O : Starting shape, n : First iteration, A : Fourth iteration, * 
; Ninth iteration. 



Figure (6) Static pressure rise through the reference diffuser at successive iterations. The grid 
size is 61 by 31. Q : Area-averaged pressure rise, A : Velocity-averaged pressure 


rise. 
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